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ABSTRACT 

This paper deals with nearly singular, possibly indefinite problems for 
which the usual multigrid solvers converge very slowly or even diverge. The 
main difficulty is related to some badly approximated smooth functions which 
correspond to eigenfunctions with nearly zero eigenvalues. A modification to 
the usual coarse-grid equations is derived, both in Correction Scheme and in 
Full Approximation Scheme. With this modification, the algorithm exhibits the 
usual multigrid efficiency. 
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INTRODUCTION 


Usual multigrid for indefinite problems is sometimes found to be very 
inefficient, A strong limitation exists on the coarsest grid to be used in 
the process. This limitation is not so much a result of the indefiniteness 
(existence of eigenvalues with different signs) itself, but of the nearness to 
singularity, that is, the existence of nearly zero eigenvalues. These 
eigenvalues are badly approximated (e.g., they may even have a different sign) 
on coarse grids, hence the corresponding eigenfunctions, which are usually 
smooth ones, cannot efficiently converge. As a remedy, one could avoid using 
grids which are too coarse, but in many cases this would degrade efficiency. 

This trouble of the coarse-grid approximation has been resolved by 
introducing a modification to the usual coarse-grid equations, based on the 
observation that there are just few smooth eigenfunctions which are not well 
represented on the coarse-grid, and these can be controlled by specially added 
relations. This modifiction removes the restriction on the coarseness of the 
grids that can be used. 

Another issue when dealing with indefinite problems is the choice of 
relaxation. Mode analysis shows that the Gauss-Seidel relaxation is suitable 
for such problems if fine enough grids are considered. Indeed, even though 
some smooth components diverge with this relaxation, on fine enough grids this 
divergence is slow and can, therefore, easily be corrected by the coarse-grid 
corrections. On coarser grids, however, the divergence of smooth components 
in Gauss-Seidel relaxation is faster, hence, another relaxation scheme is 
needed. We have used for that purpose the Kaczmarz relaxation, which always 
converges . 
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The multigrid algorithm obtained here has good asymptotic convergence 
rates for problems in which the indefiniteness is not too high, i.e., the 
number of eigenvalues with the "wrong" sign (positive in our text) is small. 
For higher indefiniteness another method has been developed and will be 
reported elsewhere. 


2. RELAXATION 

Generally, in order to achieve good multigrid performances, the relaxation 
involved need to have good smoothing properties on one hand, and at most slow 
divergence on the other hand. We discuss below Gauss-Seidel and Kaczmarz 
relaxations and their proper use in our context. 

2.1. Gauss-Seidel 

Fourier analysis of Gauss-Seidel relaxation even for slightly indefinite 

problems shows that on fine grids high frequencies converge very fast. The 

reason for this is that the principal part for indefinite problems is the same 

as that of definite ones. Smooth components may diverge on such grids, but 

slowly enough to be handled by the coarse-grid correction. For example, in 

case of the operator A + k in two dimensions the worst divergence factor 

9 9 

per sweep of smooth components is 1/(1 h ), h being the mesh size. 

On coarse grids, typically when this factor becomes larger than 1.2 or so, 
Gauss-Seidel relaxation can no longer be used. 
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2.2. Kaczmarz Relaxation 
Given an equation 

Ax = b (2.1) 

where A is an n x n nonsingular matrix, define a new unknown y such that 

A* y = x (2.2) 

where A* is the adjoint of A. The equation obtained for y is 

AA* y = b. (2.3) 

The matrix AA* is symmetric positive definite for any A which is 
nonsingular. Hence, Gauss-Seidel relaxation for equation (2.3) will 

converge. It induces a relaxation on equation (2.1) via the relation (2.2). 
This relaxation of (2.1) is called Kaczmarz relaxation. Its i-th step is 

Xj «• Xj + a^ S i (j " l,***,n) 

6 i ■ ( b i ' a lj l«y| 2 

where is the complex conjugate of a ij* For general smoothing 

properties of this relaxation, see [1] Section 1.1 and [2]. 

Kaczmarz relaxation converges whenever solution exists, and can therefore 
be used on coarse grids. Moreover, when more than one solution exists the 
convergence is to the one closest to the initial approximation (Tanabe [3]). 



Hence, this relaxation would not allow the growth of error eigenfunctions 
corresponding to X = 0, and would similarly allow only very slow change in 
eigenfunctions corresponding to X close to 0. 


3. TROUBLES WITH THE CO ARSE-GRID APPROXIMATION 

Having settled the question of relaxation, another difficulty is 

encountered: some coarse grids do not well approximate some smooth 

components. To understand this situation, suppose the error on the fine-grid 

(grid h) contains a smooth eigenfunction <j> , so that the corresponding 
h h h h 

residual is L <J> = X <(> • The corresponding equations on the coarse-grid 

(grid H) are 

l h v h = x h i“ * h 


where is the f ine-to-coarse transfer, i.e., some local averaging. Since 
h v* H Vi 

<J> is a smooth eigenfunction of L , 1^ § is approximately an 
eigenfunction of L H , but with slightly different eigenvalue \ H . The 
solution of the coarse-grid equations is approximately 


- £ i? * h - 


After interpolating the V H as a correction to the fine-grid solution, the 
new error is approximately 




*H 




"h. h 

where 1^ is the interpolation operator and since <p is a smooth function, 

we assume that ijj 1^ 4>^ « <f>^. Thus, due to the coarse-grid correction the 

error is reduced by the factor |l - X /X |; hence a condition for good 

convergence is 


|i 



(3.1) 


for any eigenfunction <)> which has poor convergence by relaxation, 
h H 

When X , X are close to zero, relation (3.1), even if it holds on fine 
enough grids, it may strongly be violated on coarse grids. Such coarse grids 
cannot then be used in the multigrid process. Without them, however, 
efficiency may very much degenerate. We will therefore present a new method 
in which restriction (3.1) is removed. 


4. MODIFIED COARSE-GRID EQUATIONS: Two-Grid Case 

The modification described here is based on the assumption that there are 
only few smooth eigenfunctions for which relation (3.1) is violated. Denote 
by Hq the subspace spanned by these badly approximated eigenfunctions. We 
assume for the description below that Hq is known. In Section 6 we present 
a method for approximating Hq. 

4.1. Correction Scheme (CS) Version 

h , 

Assume first that Hq is spanned by one function <j> , and let u h be the 
exact solution of the fine-grid (grid h) equation 


(4.1) 



^ ll 

Suppose the current approximation u to u satisfies 


/ h A hv _ v~h , .h .hv 

<U ,<()> = <U + Tpj) , (j> >• 


(4.2) 


^h h 

If r\ were known we would have the approximation u 1 + tv{> on the fine-grid 
instead of u . This would yield the coarse-grid (grid H) equation 


.H H T H r h T h ~h T h h-i 

L v = I, [F - L u - nL <(> J 


(4.3) 


H h h ~h h 

where v* 1 approximates the error v n = u - u - r\$ . Since by (4.2) the 

latter does not have components in Hq, equation (4.3) could be used to 

accelerate fine-grid convergence. However, since q is not known, we need to 

add another equation on the coarse grid which will enable us to solve also for 

n • A reasonable choice for such an equation is an approximation of equation 

(4.2), namely 

<V H , I? <j> h > = 0 (4.4) 

££ 

where 1^ is some f ine-to-coarse transfer, not necessarily identical with 
1^. Equations (4.3), (4.4) form the modified CS equations . 

Suppose now that Hq is spanned by and <<jK, = 6^* 

Because of linearity, the corrected CS equations for this case will be 


T H H t H D h 
L v = L R 
h 


l 1“ L h 
J h J 


(4.5a) 


. H H n n 
<V , <b . > = 0 


1 >N 


j - 1 


(4.5b) 
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where 



■ H _ "j® .h 

♦j " h V 


h Vi 

F - L n u , u being the current fine-grid approximation, and 
The coarse-grid correction will finally be done either by 


~h 

u 


h 

u 




~H 
v + 


N 

I 

j=i 


~ H 
T| • <p . 

J 


(4,6a) 


or by 


~h ~h 
u «- u + 


r h -H 


N 


^ + jl, "a h 


(4.6b) 


depending on whether or not <j>. are stored on the fine grid, u and p. 

J J 

being the computed (approximate) solutions to equations (4.5). The difference 

h 

between (4,6a) and (4,6b) is usually unimportant, so need not be 

stored. The only case where (4.6b) must be used is when the fine-grid problem 
is much closer to singularity than the coarse-grid one. In that case q. may 

h „ 

be large; therefore, 4>j ma Y have large high frequency components, which 

will magnify the residuals on the fine grid. By doing (4.6b) one avoids 
introducing high frequency components that arise from interpolating ?i <Jk, 
and therefore the mentioned difficulty is removed. See Section 7, Tables 5 
and 6 • 


4.2. Full Approximation Scheme (FAS) Version 

The Full Approximation Scheme is essential for nonlinear problems or when 
local refinement is used. It is important, therefore, to derive the modified 
equations in that formulation too. This derivation can be done directly from 
(4.5), but to gain an additional insight we do it independently. The usual 
FAS equation on the coarse grid is 



(4.7) 


T H H H H ,~h* 
L u = F + t, (u ) 
n 


where t^(u^) = u* 1 - lj| L* 1 u^ is the "f ine-to-coarse (defect) 

correction," F 1 * = F*\ u* 1 is the current approximation on the fine 

J.J 

grid, and 1^ is the f ine-to-coarse solution transfer (see [1, Sections 8.1- 

8 . 2 ]). 


In the present case, however, we wish to approximate on the coarse-grid 
only the part of the error which is free of Hg components; hence, the Hg 
components of the correction, £ <j\j, should be considered part of the fine 

grid approximation, replacing (4.7) by 


H H 
u 


F« + 



N 

+ I n, <|> 
j=l 3 



(4.8) 


The additional conditions, ensuring that the coarse-grid correction is indeed 
approximtely free of Hg components, can be written as 


, H H s _ ~h H s . H H. .. * 

<u , - <I h u , <j)j> + (j = 1 , • • • »N) . 


(4.9) 


Equations (4.8) and (4.9) together should determine u H and Once 
an approximate solution (u^, has been calculated, the correction 
to the fine grid solution can be done analogously to either (4.6a) or (4.6b), 
but the former option yields here a particularly simple formula, namely 


~h ~h 
u + u 


+ 


h,~H 
u 


y 


tH ~hs 

- h u ) 


(4.10a) 
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which is just the usual FAS correction formula. The latter option, which must 
be used in some extreme cases, reads 


~h ~h 
u u 


+ 



H ~ h r H W .h 

U - l T! $ J + l T\ $ . 

n j=i J J j=i J J 


(4.10b) 


Equation (4.8) is not in a form convenient for calculations on grid H. 
In case the problem is linear, a more convenient form is 


N 

H H -=H r ,H 

L u = F + I n. 

j=l J J 


(4.11) 


where 


F« = 


F + T h (u ) 


= ij R h + l h (iJ U h ) 


(4.12) 


and j ) * The solution of (4.9) and (4.11) thus involves the 2N+1 

input functions F^, \|i^, • • • <|>^> • • • »<|>^> of which should be calculated 

and stored whenever the algorithm switches from level h to level H, while 
the other 2N functions can be calculated and stored once for all. The same 
equations can be used also for nonlinear problems, but with generally 

calculated by 



r H,~h 

M u 


- .♦}) - 


H 

T h 


(u h )}. 


(4.13) 


with sufficiently small positive e. The dependence of on u is very 
crude (e.g., no dependence at all when L* 1 is linear); hence it will usually 
be unnecessary to update them on a new switch to level H. 



' - 10 - 


5. GENERAL MULTIPLE GRID EQUATIONS 

Suppose a sequence of discretization with mesh sizes h^ > h^ > ••• >h^ 
is given, where tu = 2h^ +1 . Let the h k ~grid equations be 


T k k _k 
L u = F 


where L k approximates L k+1 for k < M, and L M approximate some 
differential operator L. 

Usually, if level M well approximates the differential equation (even in 
terms of Hq components) then level M - 1 will approximate level M well 
enough for acceleration purposes. Hence, modified coarse-grid equations may 
not be needed on level M - 1. Denote by z the finest level on which 
modified equations are needed. We describe now the modified equations on 

levels k £, assuming the subspace of bad components, Hq, is spanned on 
level Z by the orthogonal set {4>j » • • • • 


5.1. CS Version 

For k < Z the equations to be solved for v k , q k on level k are 


. k k -k 


N 


L v = f - l n"l~ L— * 


k T k j £+1 x £+l 

j 


j=l j A+1 


(5.1a) 


x k ,kv k 

<v • *}> * <•! 


(j = 1 , . . . ,N) 


(5.1b) 


where 


-k T k ,<vk+l T k+1 ~k+l x 

f = W f " L v > 


(k < Z) 


(5.2a) 
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N 


<gk r ~k T k 

f = f " I tk I 
j=l J 


.£+1 £+1 
A+l L +j 


(k < A) 


(5.2b) 


j£+l _ ^A+l 


(5.2c) 


k k+1 ,~k+l A k+l s 

Pj - Pj - <V , > 


(k < A, j = 1 , • • • »N) 


(5. 2d) 


A _ n 

P J “° 


(j = 1, ••• ,N) 


(5.2e) 


«|> k = t.. <l, k+1 

k+1 


(k < £; j = !,•••, N) 


(5.2f ) 


rrk _ yk yk+1 
£+1 k+1 A+l * 


T k _ T k ,k+l 
A+l k+1 A+l 


(5.2g) 


, , i k are f ine-to-coarse grid transfers, not necessarily the same. 
k+1 k+1 


-k ~k 


k k 


, p. are the current approximation to v , p . respectively. Initial 
J J 

approximations are v k = 0, p k = 0. The input functions for level k are 

thus f k , <(> k and I k . l/ +1 <|k + 1 , (j = 1,***,N), of which only f k should be 
J J 

updated on every new switch from level k+1. 

For efficient relaxation, instead of storing f k one should store f k 

and update it whenever the p^ are changed. 

k— 1 k 

Note that p^ is designed to be a correction to p j . Thus, the coarse- 

grid corrections for 2 < k £ A will be done by the replacements 


~k ~k . ~k— 1 

n j + "j + n j 


(j = 1 > • • • »N) 


(5.3a) 


^k ^k r ~k-l T k T A+l , A+l 

f " f - ^ n j Vi L 




(5.3b) 



-12- 


~k ~k , _k ~k-l 

v «• v + I v 

k-1 


(5.3c) 


while for k = i + 1 use 


~*+l ^ ~A+1 + !*+!/-* + 


a ; ; •;) 


(5.3d) 


or 


~t + l + -Ml + jt+l ~t + \ ~l + *M > 

1 j=l J J 


(5.3e) 


(see discussion in Section 4.1 for the use of (5.3d) versus (5.3c)). 


5.2. FAS Version 


k k 


For k Z the equations to be solved for u , on level k are given 


by 


k k _ -=k r k k 
L U = F + l n ij> 

j=l J J 


(5.4a) 


z k k s k , k k , . . 

<u , <j)j> = Oj + rij cu (j N) 


where 


(5.4b) 


*) - - Vl<" k+1 > + l h *3 k+1 < *> ».5a) 


Ir If lr4- 1 1 

(hence t|)j = T k+1 ) + l^+i in t ^ ie linear case) 


^ +1 = 0 

j 


(j - N) 


(5.5b) 


k ,~k+K k rrk ~k+l _k k+1 ~k+l 

W u > - L r k+i u - Vi L u 


(5.5c) 
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. k _ A k+1 

♦d *k+l +j 


(j = 1, •**,N) 


(5.5d) 


-=k k yk ~k+l k /?;k+l k+1 ~k+lx 

F L r k+i u + W F ' L u > 


(5.5e) 


?;k -=k ^ ~k k 

F = F + I r\ Tp. 

J-l J J 


(k < £) 


(5.5f) 


?* +1 = f * +1 


(5.5g) 


k 


<j* 

^ k+1 


~k+l k s , 
u , d> .> + 
3 


^k+1 .~k+l k+l x 

a . -<u , <b . > 

3 3 


(k = A) 
(k < l) 


(5.5h) 


~k 

a. 

J 


k , ~k k 
a. + n . a. 
J J J 


(j - 1,...,N) 


(5.5i) 


k _ . k k. 
a. - <6 . , <t . > 
J J r 3 


( 3 “ 1, 


(5.5 j) 


and r\., u are the current approximations to p., u respectively. Intital 
J 3 

approximations are p_. = 0, u = Ij^ u . The input functions are F , 

k k k k ir 

<J >1 , • • • , <{>n » ^p***>^> of which only F must be recalculated each time the 

level k problem is formulated. 

“k 

For efficient relaxation, instead of storing F and o \ (j = 1,*«*,N) 

~k ~k ~k 

one should store F and o . and update them whenever ri . is updated* 

^k ^k k 

Initially (when the level k problem is set up) F = F and = ck • 

The coarse-grid corrections will be done by the replacements 


~k ~k . ~k-l 
n j * n j n j 


(2 < k < £; j = 1,...,N) (5.6a) 
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^k ^ J. V ~k-l ,k 

F + F + ) n. ib. 

L J J 


(2 < k < Z) 


(5.6b) 


~k ~k ~k-l k-1 

a . a . + n . a. 

J J J j 


(2 _< k < £; j = 1,...,N) (5.6c) 


~k ~k k r ~k-l -rk — 1 ~ks . _ . . 

U U + i C u “ ! k u J (2 < k < M). 


(5.6d) 


In case, for some 2 < k < 1+ 1, the grid k problem is much closer to 
singularity than grid k+1 problem, (5.6d) should be replaced by 


-k ~k , T k / ~k— 1 T-k-1 ~k 


N 


k-1 k-li 


-k-1 .k 


U + U + ^-l( u " *k u " l =l ) + l =l 'i 


(5.6e) 


which is the analogue of both (5,3c) and (5*3e). Of course, (5.6e) can always 
be used, but (5.6d) is somewhat simpler (cf, end of Section 4.1), 

Observe, indeed, that in the linear case 

- 4,(* k+1 ) 1- iL * k+1 = t. k $ k - I k , L t+1 ^ +1 

J k+l^j J k+1 y j £+1 

U lr rwlf+l k r. k k 

and by identifying u with I^ + j u + v + ^ the equivalence of the 


FAS and the CS is easily seen. 


5.3. Solution Process for Modified Equations 

We refer in this section to the FAS version, namely the equation 


T k k *=k , ^ k k 

L u = F + I n $ 

j-1 J J 


/ k k v k , k k 

<u , d> . > = a . + q . a . 

j j j j 


(j = 1 , • • • ,N) , 


(5.7a) 


(5.7b) 
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\c lc lc 

where the unknowns are the function u K and the constants The CS 

~k . ~k , _ 

version is treated similarly. As before, u , n.» F and a. will denote tne 

current (stored) approximations to u^, n k , the right-hand side of (5.7a) and 

the right-hand side of (5.7b), respectively. 

In relaxing equations (5.7) we distinguish between the following: 

(i) a local relaxation sweep 

Relax L k u k = F k for u k by either Gauss-Seidel or Kaczmarz, keeping 

n., and therefore also F , fixed. 

3 

(ii) a global step 

k k 

This will be the step for updating bj. , • • • ,n^> and the Hq components 
in u k by using (5.7b) together with (approximatley) the Hq components of 

A 

(5.7a). Host generally this is done by solving simultaneously for £5 ^ 

(j = 1,»«*,N) the system of 2N equations: 


<L k ju k + l ^ " <^ k + J Tij^ <1^, <j>^> (j = 1»“**N) (5.8a) 

<u k + 8, = on + q. ou (j = 1,»**,N) (5.8b) 

J J J -J V J 


and then introducing the following changes 


, . N , 

~k ~k . v ,k 

u «- u + l 8, <}>. 

j=l 3 3 


~k ~k A 
n 3 * + "j 


(j = !,•••, N) 


~k ~k y A k 

F *■ F + l n. xp . 

j=l J 3 


(5.9a) 

(5.9b) 


(5.9c) 



-16- 


~k ~k , A k 

a . = a. + n. a. 

J J J J 


(j » N), 


(5.9d) 


The local relaxation is used to smooth the error in and therefore 

should be done at all levels. On the other hand, it may be enough to do step 

(ii) on the coarsest grids only , since it deals with global variables (n^) 

and with global changes to u^. Thus, step (ii) will be done on grids k <, ra, 

where usually m < £. This will usually reduce the storage requirement of the 

k 

algorithm, since there is no need to store <(>. on levels m < k _< l. In 

k 

fact, it is often unnecessary to store even iJk for m < k < Z. Indeed for 

m < k < Z these functions are only used in the interpolation step (5.6b), 

which can be skipped in case of a V cycle, because as a smooth change to 
~k 

F , its effect on the subsequent relaxation on level k is negligible. On 
the other hand, step (5.6b) cannot be skipped in case it is followed by a 

switch back to the coarser (k - 1) grid, since in this case the smooth 

~k k 

update to F is essential . Thus, in case of W cycles, ^ must be stored 

for all levels k _< £. Generally, m < Z can be used only if no intermediate 

level k (m < k < £) is much closer to singularity than level k+1 . 


5.4. Summary. Work and Storage 

rJfc rJ\£ 

A cycle for improving u and n = (rj^ , n N ) (k O is denoted by 


r ~k ~k-> r,w-r~k ~k T k ~k ~k-, 
(u , n ) «- CMG(u , n , L , F , o J 


and is defined recursively by the following steps (A) through (D). 
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(A) Make the following v^(k) times 

(a) a local step for L k u k = F k 

(b) for k m make a global step defined in (5.8), (5.9). For k = 1, 

choose Vj(k) to guarantee convergence to small residuals, or solve the 

equations directly, and then terminate the cycle. If k > 1, continue. 

(B) Starting with u k * = * u k , * = 0 (j = 1,*»*,N) make the cycle 

K J 



k-1 

n 


«- 


CMG(u 


k-1 


~k-l 

n 



~k— 1 ~k-l \ 
F ,o ) 


y(k) times, 
k-1. 

(C) 


where F k h ck ^ are defined by (5.5) with k replaced by 


~k ~k , ~k-l 

p. *■ n. + p. 
j J J 


(k < m; j = 1, • • • ,N) 


~k ~k y ~k— 1 k 

F + F + l n. % 

j=l J J 


(k < m) 


~k ~k , ~k-l k-1 
o. . a. + p . a. 


(k < m) 


and interpolation is done either by 


~k 

u 


~k 

«- U 






or 


+ U 






k-1 


-t ' 1 


N 

- I 

J-l 


~k-l 

"j 


♦r> 


N 

+ I 

j=l 


-k-1 


V 


The second option is necessary in case the grid k problem may be almost 
singular. 
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(D) Make steps (a), (b) of (A) V 2 ^) times. 

y(k) = 2 corresponds to a W cycle; hence if y(k) = 2, the global step 
has to be done on level k, which implies m k. If y(k) = 1 we can choose 

ra to be smaller than k. Since equations (5.7) are for k _< £ < M, this 

cycle is part of a bigger cycle for k = M. 

The storage and work required by this algorithm are essentially the same 
as in the usual multigrid algorithm, since all extra work and storage involved 
are made on very coarse grids, often only on level 1, sometimes also on level 

2. In fact, £ = M-l should be used only when the finest (grid M) problem is 

itself a rather poor approximation to the differential problem, so usually 
£ < M-l, in which case the extra work is negligible compared to the work of 
relaxing grid M. 


6« APPROXIMATION OF SUBSPACE Hq 

In the preceding discussion it is assumed that Hq is accurately known. 
This section deals with how accurately Hq needs to be known and how to 
approximate it. 

6.1. Accuracy Needed for Hq 

Let ^ (i _> 1) be the smooth eigenfunctions on the finest grid 

L h <j>^ = <{>£> <<(>_£> = Sjj > (6.1) 

h h 

and let Hq be for simplicity spanned by <j>^ alone. Suppose that <f>^ is 
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not known to the algorithm, and instead <j> is used, where 


.h _ r h 

<t> “ l ^ • 

i>l 1 


( 6 . 2 ) 


Suppose an error v n = e^ <j>^ has emerged on the fine grid so that 


L h V h = e x 


(6.3) 


and the corresponding modified CS coarse-grid equations are 


T H ,.H , t H T h h . ,H 

L V + ql h L 4> = ^ ^ 


(6.4a) 


<V H , *»> = 0 


(6.4b) 


H H h H H H h 

where <}> = $ , <f>_^ = For smooth eigenfunctions <|> , we can assume 

that <j>^ are again eigenfunctions 


H H m H H 
i* 9^ ^ 9^ y 


<*l> 


(6.5) 


neglecting changes in eigenfunctions since important to our discussion are 
only changes in the eigenvalues . If we write the solution to (6.4) as 


v " I E -, ‘f’-t 


( 6 . 6 ) 


then (6.4) gives 
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X i E i + T> x i a i = 0 1 > 2 (6.7a) 

X 1 E 1 + qX i a l = X 1 e l (6.7b) 


I a E - 0. 
i>l 


(6.7c) 


Equations (6.7) imply 


n = 


q i a i e i 
6 + q x aj 


(6.8a) 


= n6/a^> = -nq i a i (i >_ 2)- 


where 


e = l q ± a il 

i>2 1 1 




x h 

1 

x H 

X 


(6.8b) 


(6.9) 


Hence, the coarse-grid correction is 


„H H 

v + n4> 


- i 


i>l 


(E ± + qa i ) ( f, E 


6 + a: 


0 + q 


1 3 1 


2 q l e l 


H , 
+ 


l 


q l a l a i 


1 ' L n 

i> 2 6 + a x 


(1 - 


. H 

q ±) e l 


( 6 . 10 ) 


H H 

Extra errors have thus been introduced in the directions of {<j> 2> <|> ...,}; 

but these should be small (relative to e^), since q A should be close to 
1 for ^ not in Hq (and also a^ will be small compared to aj by the 
condition below) and, more importantly, these errors can efficiently be 
reduced by the next coarse-grid correction. Our focus here should thus be the 
behavior of the 4>^ component. Assuming <(> E « ^ by smoothness, the 
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coarse-grid correction, when interpolated to the fine grid and subtracted from 
the old fine grid error, gives in this component the new error e^ <J>^, where 

e x = (1 - q 1 )g(e + qj a ^ ) 1 ej . (6.11) 

The main condition for convergence is therefore 




, h 2 H 

Xj a^ + Xj 6 


< 1 , 


( 6 . 12 ) 


and the convergence factor per cycle is bounded (below) by the left-hand side 

of (6.12). This bound is indeed small when X^ is a good approximation to 

Xj, i.e., when |x^ - X^| « X^. But if that is not the situation (which is 

h 

why one should want to include <j>j in Hq in the first place) then the 

ti H 

necessary condition for fast convergence is that both |x^ s| and | X^ $| 

Vi 2 

are small compared to X^ a^. Since q^ *» 1 for i 2, the condition for 
fast convergence can be summarized as 


l 

i>2 



(6.13) 


This condition implies in particular that, if X^ = 0, that is, if the 
given problem is singular, then = a^ = ••• = 0, i.e., the eigenfunction 
<j>^ must be known exactly. This seems to be too stringent, but in fact, the 
increase in accuracy for <J>^ can be obtained as the algorithm proceeds, by 
doing for each cycle of the original problem, a cycle for improving <^. 
(See Section 7, Tables 8 and 9.) 
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On the other hand, (6. 13) implies that there is no complication when the 
coarse-grid problem is singular (see Section 7, Tables 1-4). 

Generally, condition (6.13) gives a precise idea as to how closely <j>^ 
should be approximated. 


6.2 Algorithm for Approximating Hq 

We motivate the algorithm by considering the case where Hq is spanned by 
one function. It is assumed in this discussion that the finest grid problem 
is well-posed. This implies that errors in components which belong to Hq 
show sizeable average residuals on the finest grid. 

The method for approximating Hq is based on the following observation: 
components which belong to Hq are spanned by eigenfunctions whose 
eigenvalues are much closer to zero than others, and exactly such 
eigenfunctions will converge in Kaczmarz relaxation much slower than other 
eigenfunctions. Hence, if the coarse-grid equation 


L k W k = 0 


with homogeneous boundary conditions 


(6.14) 


is relaxed, starting with a random approximation, then when convergence has 
slowed down, the dominant part in the resulting W must be a component in 
Hq; therefore, W at this stage can serve as an approximation to a function 
in Hq on the coarse grid. Hq is needed on finer grids. A first candidate 
will be just an interpolation of W to these grids. However, since inter- 
polation introduces high frequency errors which will leave large residuals in 
equation (6.14) on finer grids, and therefore give the wrong ^ , one needs to 
smooth somehow the interpolated W from coarser grids. A reasonable way to 
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do it is to relax (6.14) on fine grids after obtaining a first approximation 
from coarser grids. This is summarized in the following algorithm. 

Algorithm 

Repeat the following for i = 1,»**,N (N = dim Hq) 

(A) Set k = 1 

(B) W k = random function if k = 1 

W k = I k W k_1 if k > 1. 

k-1 

(C) Relax (6.14), starting with initial approximation W , keeping 

''“’k r k k i 

W orthogonal to , • • • , } , until convergence becomes slow 

(D) k+k+1 

(E) if k £ £+1 go to (A), else = W^ + ^. 

(F) Define <f> k = I k +1 <j, k+1 (k = £, £-1 , • • • , 1 ) . 

In case N is not known in advance, stop the above procedure when step 
(C) no longer reaches slow convergence in just a few sweeps. If, after few 
cycles of solving the original problem, convergence rate still deteriorates, 
repeat (A) through (F) once more, replacing the random function in (B) by the 
residuals left by the original problem on the coarsest grid. If the addition 
of the new function to the set , ••• does not improve 

convergence rate significantly, it means that the accuracy of {<j>j + ^ , • • • , } 

is not enough and this can be improved by inverse iteration on the grid 
k = £+1 (using standard multigrid for doing the inverse iteration). The 
improvement of the functions by inverse iteration is done by one 

multigrid cycle before each multigrid cycle of the original problem. Such an 
improvement is needed when the original finest grid probelm is much closer to 
singularity than the next level; see Section 7, Tables 8 and 9. 
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Finally, if cycles on level £+1 converge satisfactorily, but not on 
finer levels, then JJ, should be increased (by 1). The above algorithm can 
then, of course, be shortened, starting with the known Hq on level £,. 


7. NUMERICAL RESULTS 

Experiments were performed with the new algorithm using the model problem 

(A + k 2 )U = F in n = (0,1) x (0,1) 

U = g on 80. 

The tables below show the residual history on the finest level. We denote 
by M, ra, 1, h^ the following: 

M — the finest level, 

SL — the finest level on which corrected equations are needed, 
m — the finest level on which the global step is performed, 
hj — the mesh size of the coarsest grid (grid 1). 

The subspace Hq was calculated by the algorithm of Section 6.2, where in 
step (B) 40 relaxations were done on the coarsest grid (k = 1) and two on 
every finer grid (1 < k £). The algorithm CMG of Section 5.3 was used with 
Vj(k) =2, V 2 (k) = 1 when Gauss-Seidel relaxation is used, 

Vj(k) = V 2 (k) = 3 when Kaczraarz relaxation is used, 
v x (k) = 13 for k = 1, 

y(k) = 1 for k > 2, y (k) = 2 for k = 2. 
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We give below the discrete first two eigenvalues of the Laplacian for the 
grids used in the examples of this section. This will enable us to see how 
far from singularity is each of the different levels used in the process. 


h 

x h 

A 1 

A 2 

.25 

-18.74516600406 

-41.37258300203 

.125 

-19.48683967711 

-47.23375184668 

.0625 

-19.67587286709 

-48.81161578777 

.03125 

-19.72335955067 

-49.21342550952 


In Tables 1-4 interpolation of corrections was made according to (5.6b); 
in other tabels the interpolation is specified. Residuals were transferred by 
9-point full weighting and the local relaxation was Gauss-Seidel for kh < .5 
and Kaczmarz for kh > .5. In all examples M = 4, h^ = .25, £ = 3, m = 2. 

Tables 3, 4 show a case in which the second eigenvalue is very close to 
zero, and its corresponding eigenspace is two-dimensional. Therefore, only 
two functions were used in spanning Hq. The algorithm for finding functions 
in Hq finds first these eigenfunctions whose eigenvalues are closest to 
zero. Therefore, the eigenfunction belonging to the first eigenvalue was not 
used in these computations, and it was not needed as can be seen from the fast 
convergence shown by these tables. This clearly shows that Hq is related to 
almost-singularity , not to indefiniteness. 

In Tables 5, 6 we show that in case the finest grid problem is too close 
to a singularity one must use interpolation of correction according to (5.6e) 
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(or (5.3e) in the CS version) and not the usual FAS interpolation. In these 
two tables the exact Hg was used. 

Table 7 shows that if <f> is not known accurately enough, poor results are 
obtained. 4> in this example is found by the same procedure described in the 
beginning of this section. 

In Table 8 inverse iteration (done by usual multigrid) was used to improve 
the accuracy of <)>. Starting with $ obtained as in Table 7, one multigrid 
cycle of inverse iteration was done to improve cf> before each multigrid cycle 
for the original problem* Results are identical to the ones obtained with the 
exact <(> (Table 5). 

Table 9 shows a case in which the distance of the closest eigenvalue to 

—8 

zero is about 1.10 . As seen from this table, improving <j> by only one 
cycle of inverse iteration per cycle of the original problem is not quite 
enough to maintain the full speed of the algorithm. Once in few cycles the 
residuals are magnified, and this happens whenever the L. 2 ~norm of the error in 
the approximation is reduced significantly. This reduction of the error is 
due to a correction of the approximate solution by r\c[). If <{> is not 
accurate enough, components other than the desired ones enter to u^ and 
since their residuals are much higher than those of a magnification of 
the residuals occurs. (A similar phenomenon can be seen also in Table 7, 
where the distance from singularity is larger, but <}> is not improved at all 
by inverse iterations.) This would not have happened if we allowed the speed 
of convergence of the inverse-iteration cycles to be slightly faster than that 
of the main cycles, e.g., by adding an extra inverse-iteration cycle once per 
several cycles. But this is not really needed, because all that may happen is 
a minor slowdown at high-accuracy solutions (much below truncation errors) for 
cases of extreme closeness to singularity. 
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In fact, we believe that, to obtain solutions with errors smaller than 
trunction error, all one has generally to do is a one-cycle FMG algorithm for 
calculating Hq (meaning one inverse-iteration cycle for level k after step 
(C) in the algorithm of Section 6.2), followed by a one-cycle FMG algorithm 
for solving the original problem. 
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Tatele i ; k 2 = 18.745166, dim Hq = 1 


cycle // 

II residuals II ^ 

i 

.363 (+3 ) 

2 

,172(+2) 

3 

•114(+1) 

4 

.891(-1) 

5 

. 7 62 (—2 ) 

6 

.685(-3) 

7 

.652(-4) 

8 

.658C-5) 

9 

,684(-6) 

10 

• 7 44 (—7 ) 


Table 2: k 2 = 

19.486839, dim H Q = 1 

cycle # 

II residuals II ^ 

i 

,363(+3) 

2 

. 174(+2) 

3 

.114(+1) 

4 

,892(-l ) 

5 

. 763(-2) 

6 

.687 (-3) 

7 

.654(-4) 

8 

.661 (—5 ) 

9 

,688(-6) 

10 

• 7 49 (—7 ) 
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Table 3: k 2 = 41.372583, dim H Q = 2 


cycle # 

II residuals II 2 

1 

,363(+3) 

2 

,l72(+2) 

3 

. 112C+1) 

4 

.938(-l) 

5 

,864(-2) 

6 

,832(-3) 

7 

. 820 (— 4 ) 

8 

,815(-5) 

9 

.796(-6) 

10 

•811 (—7) 


Table 4: k 2 = 

47.233752, dim H Q = 2 

cycle # 

II residuals II 2 

1 

.363(+3) 

2 

.171(+2) 

3 

.110(+1) 

4 

• 9 1 0 (—1 ) 

5 

,824(-2) 

6 

.778(-3) 

7 

.755(-4) 

8 

.740(-5) 

9 

,682(-6) 

10 

.673(— 7) 
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Table §: 


Table 6: 


- 19.723368, dim Hq = 1, interpolation according to (5.6e). 


2 = 


cycle # 

II residuals II ^ 

1 

,363(+3) 

2 

. 172(+2) 

3 

. 114(+1 ) 

4 

• 893 (— 1 ) 

5 

,765(-2) 

6 

,687(-3) 

7 

.691 (—4 ) 

8 

.685(-5) 

9 

.759(-6) 

10 

■ 

00 

/^s 

1 

19.723368, dim Hq = 

1, interpolation acc 

cycle # 

II residuals II ^ 

1 

.363 (+3 ) 

2 

. 172(+2) 

3 

.114 (+1 ) 

4 

.893 (— 1 ) 

5 

.125 

6 

. 102 (—1 ) 

7 

. 162(-1 ) 

8 

. 132C-2) 

9 

.909 (—3 ) 

10 

.742 (—4 ) 
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Tafele 7: 


Table 8: 


k.^ - 19.72336843, dim Hq = 1, interpolation according to (5.6e), 
(j)^ crudely computed. 


cycle // II residuals II 2 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 


,363(+3) 
. 174(+2) 
.114(+1) 
.879 ( — 1 ) 
.116 

. 5 50 ( —2 ) 
.134 

.427 (-1 ) 
. 465 (— 2 ) 
.534(-3) 


= 19.72336843, dim Hq = 1, interpolation according to (5.6e). 
<j>j successively improved. 


cycle # II residuals II 2 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 


.363 (+3 ) 
. 174(+2) 
.114(+1) 
,893(-l) 
.765(-2) 
.687 (-3) 
.691(-4) 
.683(-5) 
♦759(-6) 
.768(-7) 
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Tahl@ 9: = 19.72335955955, dim Hq = 1, interpolation according to (5.6e). 


cycle // 

II Residuals 11^ 

;~h h„ 

llu - u II 2 

i 

»363(+3) 

.555 

2 

.174(+2) 

.392 

3 

•114(+1) 

.392 

4 

.893(-l) 

.392 

5 

.764(-2) 

.392 

6 

.687(-3) 

.392 

7 

.655(-4) 

.392 

8 

.359(-3) 

.268(-l) 

9 

•284(-4) 

.268(-l) 

10 

,219(-5) 

,268(-l) 
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